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Abstract 

We study the interplay of hydrodynamic mesoscale structures and 
the growth of plankton in the wake of an island, and its interaction with 
a coastal upwelling. Our focus is on a mechanism for the emergence 
of localized plankton blooms in vortices. Using a coupled system of a 
kinematic flow mimicking the mesoscale structures behind the island 
and a simple three component model for the marine ecosystem, we 
show that the long residence times of nutrients and plankton in the 
vicinity of the island and the confinement of plankton within vortices 
are key factors for the appearance of localized plankton blooms. 
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1 Introduction 



The interplay between hydrodynamic motion and the distribution of marine 
organisms Uke phytoplankton and zooplankton is a ma j or challenge recently 



addressed in numerous studies (IMann and Lazier .ligOlMDenman and Gargett 



19951: Abrahaml.ll9 98: Pet ers and Marrasel.l200(]l:lKarolvi et al 



2000; 



2001al lbl: iMartin et al., . ,20021 : iMartinl . l2003l : ISandulescu et al.l . l2007bi ^ . 



Lopez et al 



The growth of phytoplankton in the world's oceans depends strongly on 
the availability of nutrients. Thus, one of the essential factors controlling 
the primary production is the vertical transport of nutrients. Coastal up- 
welling is one of the most important mechanisms of this type. It usually 
occurs when wind-driven currents, in combination with the Coriolis force, 
produces Ekman transport, by which surface waters are driven away from 
the coast and are replaced by nutrient-rich deep waters. Due to this nutrient 
enrichment, primary production in these areas is strongly boosted, giving 
rise also to an increase of zooplankton and fish populations. 

On the other side, the interplay between plankton dynamics and hori- 
zontal trans port, mixing and stirring has been investigated in several stud 
ies recentlv (lAbrahaml. Il998l: 



2002: Hernandez-Garcia et al 



Lopez et al.l. l2001bl: 



2003; 



Hernandez-Garcia et al 



Martini . |2003| ). Horizontal stirring by 



mesoscale structures like vortices and jets redistributes plankton and nutri- 

ents and may enhance primary production ( Martin et al. . 20021 ; Hernandez-Garcia and Lopez , 
200i). Horizontal transport can also initiate phytoplankton b looms and af- 



fects competition and coe xistence of different plankton species (jKarolvi et al. 
2000l : iBracco et all , bood ). 



Vertical upwelling in connection with strong mesoscale activity occurs 
in several places on Earth. One of these regions is the Atlantic ocean area 
close to the northwestern African coast, near the Canary archipelago. The 
main water current in this area flows from the Northeast towards the Ca- 
nary isl ands, in which wake s trong mesoscale hydrodynamic activity is ob- 



served ( Aristegui et al.l . ll997l l. The interaction between the vortices emerg- 



ing in the wake of the Canary islands and the Ekman flow seems to be 



jroduction in the ope n 
20041 ). 



essential for the observed enhancement of biological 
southern Atlantic ocean close to the Canary islands (jAristegui et al. 
The aim of this paper is to study the interplay between the redistribution 
of plankton by the vortices and the primary production. In particular we 
focus on the role of residence times of plankton particles in the wake of the 
island. Though we believe that our study is relevant for different areas in 
the world, we focus on the situation around the Canary archipelago to be 
specific. 
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In this work we consid er the coupHng of the kinematic flow introduced in 
( Sandulescu et al. . 2006al ) to a simphfied model of plankton dynamics with 
three trophic levels, and study the impact of the underlying hydrodynamic 
activity and the upwelling of nutrients on primary production in different 
areas of the wake. In this setup vortices have been re ported to play an 



essent ial role in the enhancement of primary production (j Sandulescu et al 



Our main objective here is to analyze this mechanism in detail 



and show that the extended residence times of plankton within vortices are 
responsible for the observation of localized algal blooms in them. 

The organization of the paper is as follows. In section [2] we present the 
general framework of our system, indicating the hydro dynamical and the 
biological model, as well as their coupling. Our main analysis is devoted 
to the mechanism of the appearance of a localized plankton bloom within a 
vortex (Sec. [3]). We study the residence times of plankton within vortices 
and in the neighborhood of the island. Additionally we clarify the role of 
the chaotic saddle embedded in the flow in the wake of the island. Finally 
we summarize and discuss our results in Sec. IH 



2 General framework: velocity field, plankton model 
and boundary conditions 

Our system consists of a hydrodynamic flow with an embedded obstacle and 
vortices in its wake. The model contains also a current perpendicular to the 
main flow that models an Ekman flow coming from the coast, and a nutrient- 
rich region at a distance from the obstacle simulating a coastal upwelling 
zone. A sketch of the model is shown in Fig. [H With this simplified geome- 
try we mimic the essential features of the hydrodynamic flow in the Canaries 
(note that the whole Canary archipelago is approximated by one cylindrical 
island). In particular, in the wake of the obstacle strong mesoscale activity 
is observed in the form of a periodic detachment of vortices, which then 
travel in the main flow direction. 



We use the kinematic model first developed bv lJung et al.l (119931), which 
we modified by the introduction of the Ekman flow ( Sandulescu et al. . 2006al ). 
This model is coupled to a simple population dynamics which features the 
interaction of nutrients N, phytoplankton P and zooplankton Z. The next 
two subsections are devoted to the introduction of the hydrodynamic as well 
as the biological model before discussing the results of coupling both models 
to study the feedback between hydrodynamics and phytoplankton growth. 
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Figure 1: The simplified island wake model setup. 



2.1 The hydrodynamic model 



We no w introduce the velocity field. Details can be found in (jSandulescu et al.l . 



2006al l. The setup of our hydrodynamic model is based on a horizontal flow 



pattern. As Fig. [T] shows the main current runs from left to right along the 
horizontal x direction. The center of the cylinder is placed at the origin of 
the coordinate system. We consider a two-dimensional velocity field which 
can be computed analytically from a stream function The velocity com- 
ponents in X- and y-direction and the equations of motion of fluid elements 
are: 

d 

X = Vx{x,y,t) = —'i?{x,y,t), 
oy 

d 

y = Vy{x,y,t) = -—'^{x,y,t). (1) 



The stream function is given by the product of two terms (jJung et al.l . 

19931 ): 
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^{x,y,t) = f{x,y)g{x,y,t). (2) 

)undary conditions at the 
The second factor g{x, y, t) 



The first factor f{x,y) ensures the correct boundary conditions at the 

, 2" 



cyhnder, /(x, y) = 1— exp —a yyx^^+y^ — 1 

models the background flow, the vortices in the wake, and the Ekman flow 
g{x,y) = -whi(t)gi{x,y,t) + wh2{t)g2{x,y,t) + uos{x,y)y + ueQ{x - l)x. 
The vortices in the wake are of opposite sign but their maximal vortex 
strengths are equal an d denoted by w, and its sh ape is described by the func- 
tions gi (see details in Sandulescu et al. ( 2006al )). The characteristic linear 

— 1/2 

size of the vortices is given by Kq and the characteristic ratio between the 
elongation of the vortices in the x and y direction is given by a. The vortex 
centers move along the x direction according to xi{t) = 1 + L (t/Tc mod 1) 
and X2{t) = xi{t — Tc/2), and at values ofyi{t) described below. Each vortex 
travels along the x direction for a time Tc and disappears. The background 
flow moves in the positive horizontal direction with a speed uq. The factor 
s{x,y) describes the shielding of the background flow by the cylinder in a 
phenomenological manner, using the same elongation factor a as for the 
vortices. The Ekman drift, which is intended to model the flow from the 
coast towards the ocean interior, is introduced by considering an additional 
velocity of constant strength w^; in the y direction acting only at x coordi- 
nates larger than 1, i.e. just behind the island. This corresponds to a stream 
crossing the vortex street towards negative y values beyond the cylinder. 

Real oceanic flows are never perfectly periodic. Therefore we use a non- 
periodic version of the kinematic flow just presented. Non-periodicity is 
achieved by adding some randomness to the vortex trajectories. Instead of 
moving along straight horizontal lines, yi{t) = yo, y2{t) = —yo {yo constant), 
the vertical coordinates of the vortices move according to yi{t) = yo + ^Ci't^)^ 
and 2/2 (i) = —yiit), where, at each time, ^(i) is a uniform random number 
in the range [—1,1], and e is the noise strength. 

The parameters of the model are chosen in such a way that they rep- 
resent properly the geophysical features of the Canary zone. These values 
are given in Table 12. 1[ To make the model dimensionless we measure all 
lengths in units of the island radius r = 25 km and all times in units of 
the period = 30 days . For a discussion of all parameters we refer to 
( Sandulescu et al. . 2006al ). where the adaptation of the model to the situa- 
tion around the Canary islands is discussed in detail. 
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Table 1: List of parameters used in the hydrodynamical model 



parameter 


value 


dimcnsionlGSS value 


r 


25 km 


1 


Uo 


0.18 m/s 


18.66 




25 km 


1 


a 


1 


1 


w 


« 55 X 10^ Tin? /s 


200 


Tc 


30 days 


1 


L 


6r = 150 km 


6 


a-i/2 


25 km 


1 


UE 


0.02 m/s 


2 


yo 


r/2=12.5 km 


0.5 


e 


6.25 km 


0.5 



2.2 The biological model 

One can find in the literature a large variety of different models used to anal- 
yse the dynamics of marine ecosystems. Their c omplexity ranges from s imple 



ones with only a few interacting co mponents (ISteele arid Hendersonl . Il981 



19921 ) to large comprehensive ones (jBaretta et al.l . Il997 ). We use a system 



Steele and Henderson 



which is based on a three com ponent model developed by 

(119921 ) and later modified by ( Edwards and Brindleyl . ll996l ) and (|Oschlies and Garcon 
I1999I ). 

The model describes the interaction of three species in a trophic chain, 
namely nutrients A^, phytoplankton P and zooplankton Z, whose concen- 
trations evolve in time according to the following dynamics: 



dN 



dP_ 

m 

dZ_ 
Itt 




a + rjP^ 



(3) 



Le t us n ow briefly discus s the m eaning of the different terms (cf . iQschlies and Gargon 
( I999I ) and Pasquero et al. ( 2004 ) for details): the dynamics of the nutri- 
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ents is determined by nutrient supply due to vertical mixing, recycling by 
bacteria and consumption by phytoplankton. Vertical mixing which brings 
nutrients from lower layers of the ocean into the mixed layer is parameter- 
ized in the biological model, since the hydrodynamical part considers only 
horizontal transport of nutrients. For the vertical mixing we assume Nq as 
a constant nutrient concentration below the mixed layer. Thus the mixing 
term reads: 



^N = S{x,y){No-N), (4) 

where the function S determines the strength of the upwelling and will 
be discussed in more detail below. The nutrients are consumed by phyto- 
plankton with a saturation characteristics described by a HoUing type II 
functional response. The recycling by bacteria is modelled by the last three 
terms in the bracket. A part of all dead organic matter as well as the exu- 
dation of zooplankton is degraded by bacteria, though the dynamics of the 
bacteria themselves is not included in the model. The phytoplankton grows 
upon the uptake of nutrients, but its concentration is diminished by zoo- 
plankton (grazing term) and due to natural mortality. Grazing, modelled 
by a Holling type III function, enters also as a growth term for the zoo- 
plankton dynamics multiplied by a factor 7 taking into account that only 
a part of the food is converted into biomass of the zooplankton, while the 
other part (1 — 7) goes to recycling. The natural mortality of zooplankton 
is assumed to be quadratic because this term does not only model natural 
mortality b ut also the existence of hi gher predators which are not explicitly 
considered ( Edwards and Beed . 200ll ). The parameters used are taken from 



(jPasquero et al. . .2004 ) as presented in Table 12.21 Although appropriate for 
open ocean, they would provide estimates for biological properties in the 
Atlantic not too close to the coast. To obtain dimensionless quantities con- 
venient for the numerics, space is measured in units of r, time in units of Tc 
and mass in units of lO^'^ mmol N. 

The primary production is defined as the growth term in the phytoplank- 
ton dynamics: 

The function S, measuring the strength of vertical mixing in this model 
is a crucial quantity for the coupling between the hydrodynamical and the 
biological model, because it quantifies the local nutrient supply. As shown 
in Fig. [T] we assume that there exists an upwelling zone which is located 
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Table 2: List of parameters used in the biological model 



parameter 


value 


dimensionless value 


/? 


0.66 day-i 


19.8 


V 


1.0 (mmol N m^"^)^^ day^^ 


0.12288 


7 


0.75 


0.75 


a 


2.0 day-i 


60 


Si 


0.00648 day~^ (nutrient poor) 


0.1944 


Sh 


0.648 day~^ (nutrient rich) 


19.44 


kN 


0.5 mmol N m~^ 


7.8125 


^J^■N 


0.2 


0.2 


HP 


0.03 day 


0.9 




0.2 (mmol N m~^)^^ day^^ 


0.384 


No 


8.0 mmol N m~^ 


125 



in a small rectangular region on one side of the island. According to this 
assumption, we assign two different values to the parameter S{x,y). In 
the upwelling zone there is a strong vertical mixing leading to nutrient rich 
waters in the mixed layer. There we assume S{x,y) = Sh = 0.648 day""*^, 
while in all the surrounding waters upwelling is much lower so that we assign 
S{x, y) = Si = 0.00648 day-^ which is a hundred times lower. 

The dynamics of this model is different depending on the choice of 
parameters. The long-term behavior can be either stationary with con- 
stant c oncentrations of N, P a nd Z or os cillatory. We re f er for more de- 
tails to Edwards and Brindley ( 19961 ) and Pasquero et al. I (120041 ). We use 



a parameter set where the system possesses a stable steady-state. Us- 
ing the parameter values from Table 12.21 and setting the vertical mixing 
S = Si = 0.00648 day~^ we obtain as a steady state Namb = 0.185, 
Pamb = 0.355 and Z = Zamb = 0.444 mmol N m~'^. In this nutrient poor 
region the ambient primary production is PPamb = 0.0633 mmol N m~^ 
day-^ 

2.3 The coupled model 

The coupling of the biology and the hydrodynamics yields a system of 
advection-reaction-diffusion equations. Thus the complete model is given 
by the following system of partial differential equations: 
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dN , 

+ V.S/N = Fn + DV^N, 

ot 
dP 

'dt 



+ V-VP = Fp + DV^P, 



dZ 

— +V-VZ = Fz + DV^Z, (6) 

with the biological interactions F/v, Pp, and Fz from Eq. and the 
velocity field v{x,y,t) from Eqs. ([T|) and The diffusion terms take 

into account the small scale turbulence with eddy diffusivity D. We take 
D ~ lOm^/s, as correspo nding to the Okubo estimation of eddy diffusivity 



at scales of about 10 km (|Okubd . Il97lh . the scales which begin to be missed 



in our large scale streamfunction. This advection-reaction-diffusion system 
is solved numerically by means of the method explained in Appendix A. As 
we are studying an open flow, the inflow conditions into the left part of the 
computational domain have to be specified to complete the model definition. 
Depending on the choice of the inflow concentr a,tions we observe differen t 
behavior. A detailed analysis can be found in ( Sandulescu et al. . 2007bl ) . 



here we only recall the main results which are the basis of the analysis we 
present here. We have studied two essentially different inflow conditions: 

1. In the first one fluid parcels enter the computational domain with the 
ambient concentrations Namh-, Pamb, and Zamb corresponding to the 
steady-state for S = Si. In this case the exterior of the computational 
domain has the same properties as the part of the domain without 
upwelling. 

2. In the second one fluid parcels transported by the main flow enter 
the domain from the left with a very small content of nutrients and 
plankton, corresponding to a biologically very poor open ocean outside 
the considered domain. In particular we take Nl = 0.01 x Nambi 
Pl = 0.01 X Pamb and Zl = 0.01 x Zamb-, leading to very low primary 
production in the inflow water PPl = 8.6 x 10~^ mmol N m^^ day^^. 
Since those concentrations are very low, we take into account that 
fluctuations may be important by adding to each of the concentrations 
{Nl, Pl, Zl) some Gaussian noise term (cf. Sandulescu et al. ( 2007bl ) 
for details). 

The two inflow concentrations yield different behavior as shown in Fig. 
[2j In the first case we observe high biological activity connected with a high 
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primary production in the area outside vortices. Namely, this is the area of 
the nutrient plumes advected from the upwelling region (Fig. [2] left column). 
By contrast, in the second inflow case (Fig. [2] right column) we obtain a 
high phytoplankton concentration within the vortices. It turns out that here 
the vortices act as incubators for primary production leading to localized 
plankton blooms. In the first inflow case the behavior is easy to understand 
since due to higher nutrient concentrations in the upwelling region and its 
neighborhood a high growth of phytoplankton is expected. The response on 
the upwelling of nutrients in the second case is less obvious. Therefore the 
main objective of this paper is to find out the mechanism of the localized 
plankton blooms. In the rest of this work we consider only the situation 
that the concentrations at inflow are at their low values Nl,Pl, Z^. 



3 The mechanism of emergence of locahzed plank- 
ton blooms 

After specifying the complete model and its dynamics we now investigate 
the behavior of the coupled biological and hydrodynamical system from 
different perspectives to clarify the mechanism of localized enhancement 
of phytoplankton and primary production connected to vortices. Firstly 
we study the biological model alone to understand the interplay between 
the three biological components A'^, P and Z leading to a sharp increase of 
phytonplankton for some time interval. This study yields a certain biological 
time scale for the growth of plankton which we compare in a second step to 
the hydrodyamical time scale obtained from the investigations of residence 
times in vortices. Thirdly we discuss the role of the chaotic saddle embedded 
in the flow for the emergence of localized enhanced plankton growth. 



3.1 Plankton growth 

To study the enhancement of primary production and the emergence of lo- 
calized algal blooms we have first to analyze the dynamics of the biological 
model. There is no commonly accepted definition of an algal bloom. Usu- 
ally a large increase in the phytoplankton concentrations is considered as a 
bloom. In most cases such blooms are observed once or twice a year due to 
seasonal forcing. In our case the phytoplankton bloom is not related to an 
external forcing and appears only on a rather short time scale. We consider 
the case where there appears a sharp increase in phytoplankton as a result of 



an enrichment with nutrients (jEdwards and Brindlevl . Il996l : iHuppert et al.l . 
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20021 ). 

Since the long-term behavior of our model is stationary for the parametriza- 
tion used, the emergence of a sharp increase in phytoplankton is a transient 
phenomenon and its time scale is important for the mechanism of localized 
enhancement of the primary production. The time evolution of the three 
components and the primary production PP of the model system towards 
the steady state concentrations Camb (C = ^, P, Z, or PP) is shown in 
the upper part of Fig. [3l With starting concentrations 0.01 x steady-state 
concentrations of Namb, Pamb and Zamb, first the nutrient concentration 
increases and, after a time lag, primary production and phytoplankton con- 
centration follow with a large increase. This growth is approximately expo- 
nential when the nutrients reach their maximum. Finally, with a larger time 
lag the concentration of predators (zooplankton) increases as well and the 
bloom ends due to two factors: nutrient depletion and increased grazing by 
zooplankton. For comparison in the lower panel of Fig. [3] the time evolution 
of the system with starting concentrations 0.3 x steady-state concentrations 
of Nambi Pamb and Zamb is plotted. With higher starting concentrations 
the overshooting in nutrient and phytoplankton concentrations at the be- 
ginning of the time evolution is less pronounced (because there are more 
predators already present) and the concentrations converge faster towards 
the steady-state. 

From these simulations we can estimate the time scale for the biological 
growth: To reach the maximum of the bloom, only 15 to 25 days are neces- 
sary depending on the initial condition. The time scale for the whole relax- 
ation process is about 2Tc, i.e. about 60 days. To understand the interplay 
between the biological growth and the hydrodynamic mesoscale structures 
we have to compare this biological time scale with the hydrodynamic one. 



3.2 The residence time of fluid parcels in the wake 

As pointed out in ISandulescu et al.l (j2007bl ) the hydrodynamic mesoscale 
structures are important for the enhancement of primary production in the 
wake of the island. To gain more insight into the interplay of hydrodynam- 
ics and plankton growth we now quantify the time scales for the relevant 
hydrodynamic processes. To this end we study the various structures in 
the hydrodynamic flow which have a significant influence on the residence 
times of nutrients and plankton in the wake of the island. Firstly, far away 
from the island (top and bottom of Fig. [T]) the flow is strain dominated and 
particles like nutrients and plankton are advected with the background flow 
of speed uq. Thus the residence time of particles released away from the 
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island (with y > 2 and y < —2, a; = 0) is about 16 days. 

Secondly we note the existence of the eddies. They are characterized by 
a dominance of vorticity compared to strain. Thus particles are trapped in 
the vortex once entrained to it. The particles will rotate in the vortex for 
some time, but since this confinement is not perfect and vortices exist only 
for some time they leave the vortex and move away with the background 
flow out of the computational area (cf. Fig. H]). 

Thirdly we consider two other geometrical objects which are also relevant 
for the residence time of particles in t he vortex street: the chaotic saddle and 
the c ylinder boundary. As shown in (|jung et al.1 . EJ : IPuan and Wiggind . 



19971 ) there exists a chaotic saddle which is embedded in the flow beyond 



the island. At least for short time scales, this invariant set determines the 
residence times of particles. Particles released in the neighborhood of the 
chaotic saddle will approach it along its stable manifold and will leave it 
along its unstable manifold. The unstable manifolds at two different times 
(the manifolds and the saddle are tim e-dependent ) are p resented in Fig. [5l 
As compared with the configuration in Jung et aP ( 19931 ). the manifolds are 
tightly packed close to the cylinder, because of the parameters used here. 
To obtain an estimate for the residence t ime on the chaotic saddle we use 
a method suggested by Jung et aP (jl993l ). We sprinkle a large number N 
of tracer particles (A^ = 62500) in the area [0,2] x [—1,1] and integrate 
their trajectories forward in time. If the dynamics in the region is mainly 
hyperbolic, the number of particles remaining in the area of the chaotic 
saddle decreases as A^ ~ exp(— k t) with k, the escape rate or 1/k the mean 
residence time on the saddle. Figure [6] shows the residence times obtained 
with this method, and Fig. [7] shows the decay of the number of particles in 
the region as a function time. We note that the expected exponential decay 
occurs only for very short time scales. By fitting this initial time decay, the 
corresponding escape rate is k = S.l/Tc, and therefore the residence time of 
tracers in the hyperbolic part of the saddle is r 10 days. For larger times 
the particle number in the region decays as a power law. The reason for this 
power-law behavior is the non-hype rbolic dynarnics n ear the boundary of 
the cylinder. As already shown by ( Jung et al. . 19931 ) particles stay for a 
long time in the vicinity of the cylinder giving rise to another long-term 
statistics of the residence times of the tracers. Thus the number of particles 
decays as A^ ~ t~'^ with 7 = 0.96. The residence times in the vicinity of 
the island can be estimated as r ~ 85 days, measured from the decay to a 
fraction e of the initial number: N{t) = N(0)/e. 

Overall we obtain a residence time statistics which reflects a combination 
of the three components in the flow: the cylinder, the chaotic saddle and 
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the vortices. This leads to residence times of tracers in the wake as long as 
90 days (cf. Fig. [6]). Note that in Fig. [6] one can see that tracers having 
the longer residence times are either located close to the cylinder or on the 
chaotic saddle. The residence times in the vortices are determined by their 
travel time which is about 50 days. 

3.3 The interplay of biological and hydrodynamical residence 
times 

To understand the emergence of localized enhancement of primary produc- 
tion we have to analyze the interplay of the different time scales relevant 
for coupled biological and hydrodynamical processes. Biological evolution 
needs about 30-60 days to reach the steady-state when entering the com- 
putational domain with very low concentrations of nutrients and plankton. 
Due to the exponential growth in the beginning of the growth phase, we 
obtain a plankton bloom after about 25 days. Outside the vortex street the 
travel time of tracers through the computational domain is only about 16 
days due to the background flow of uq = 0.18 m/s. Therefore we cannot 
expect a considerable growth of plankton outside the vortex street, since the 
residence time of plankton and nutrients is too short. 

Let us now analyze the situation within the mesoscale structures of the 
flow in the wake of the island. As the residence time close to the island 
is about 85 days the concentrations of nutrients and plankton have already 
reached the steady-state which is also indicated by the green color in the 
right column of Fig. [2j Some of the particles in the vicinity of the island 
come close to the stable manifold of the chaotic saddle visible as the filaments 
which detach from the cylinder. These filaments are stretched and folded 
along the unstable manifold of the chaotic saddle, being diffusively diluted 
during the process by mixing with the poor surrounding waters. Thus, 
very thin filaments of low plankton and nutrient concentration are produced 
which are first rolled around the vortices and then entrained by them. Inside 
the vortex the concentrations become homogeneised to a low value. These 
very low concentrations of plankton experience the bloom cycle described 
in Section [3.11 during the time they are trapped and confined by the vortex. 
Since the travel time for the vortices is about 50 days, plankton in them 
has time to grow. Therefore we observe a localized plankton bloom when 
the vortex has travelled a distance of ~ 100 km which corresponds to a 
residence of the plankton in the vortices of ~ 15 — 20 days. After 40-50 
days and ~ 200 km we obtain steady-state concentrations and the former 
filamental structure within the vortex is smeared out by our diffusion term 
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which mimics smah scale turbulence. 



3.4 The emergence of filamental structures due to strong 
mixing 

In the previous subsection we have stated that the transport of nutrients and 
plankton from the vicinity into the interior of a vortex happens by filaments 
which are entrained by the vortex. To explain this stretching mechanism 
we now study the mixing process around the vortices in more detail using 
a method to visualize exponential divergence of the trajectories of initially 
nearby particles. 

The usual tool to analyse exponential divergence in dynamical systems 
theory is the computation of Lyapunov exponents. In order to adjust this 
concept to local processes, we compute finite size Lyapunov exponents (FSLE) 
which are based on the idea that one measures the time necessa ry to obtain a 
final prescribed distance 6f starting from an initial distance 6q ( Artale et al.l . 



19971 : Id'Ovidio et aP . hooi ) . For a two-dimensional flow we obtain two Lya- 



punov exponents A+ and A_ (see appendix B). 

Maxima in the spatial distribution of A+, the positive or exp anding 



FSLE, approximate the und erlying stable manifold of the chaotic flow (jJoseph and Legras 



2002; Id'Ovidio et al.l . l2004l ) , the direction along which parcels approach the 
saddle. The contracting FSLE, A_, detects the underlying unstable manifold 
of the chaotic flow, the direction along which parcels are stretched out of 
the saddle. For details on how to calculate these scalar fields see Appendix 
B. 

The FSLEs were calculated choosing the initial separation Jq equal to 
the gridsize and the final separation 5f equal to the radius of the island and 
the vortices, since this is the scale of the motion in the wake. As both stable 
and unstable i nanifolds cannot be c r ossed by fluid parcels they are barriers 
(|Artale et all . Il997l : Id'Ovidio et ID . l2004l ). The scalar field A+ - |A_| is 



plotted in Fig. [HI FSLEs are Lagrangian measures, which are computed 
from trajectories that remain in the fiow for a long time, in our case for up 
to 3Tc. Therefore even though they are plotted as a snapshot, the visualized 
structures reflect the stretching and folding of the fluid parcels during this 
long time. The stable and unstable manifolds are intertwined around the 
vortex cores and at the island. Stable and unstable manifolds are crossing 
the wake allowing for transport across the vortex street. They intersect 
each other in hyperbolic points, regions of strong mixing. This stretching- 
compressing mechanism leads to low nutrient and plankton concentrations 
transported into the interior of the vortex, and thus becoming the starting 
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concentrations for the focalized plankton bloom. 

3.5 On the role of the upwelhng region of nutrients 

Finally we discuss the importance of the vertical mixing of nutrients in the 
upwelling zone for the emergence of a plankton bloom inside vortices. Com- 
paring Fig. [2] left and right column it is obvious that in the case of an 
inflow with steady-state conditions (left column), the nutrient plume which 
appears in the neighborhood of the upwelling zone gives rise to a phyto- 
plankton bloom (red filamental plume). Such a plume is almost absent 
under low inflow conditions (right column). Though the nutrient supply 
due to vertical mixing is identical for both inflow conditions, it seems to 
have a limited effect in the low inflow case. One argument has been already 
discussed above: The background flow transports the nutrients too fast so 
that the very small plankton concentrations can not grow to reach high val- 
ues during the travel time through the computational area. The growth of 
phytoplankton is visible only further downstream. This leads to the con- 
clusion that the plankton bloom inside the vortex is only slightly influenced 
by the extra nutrients entrained from the upwelling zone in the low inflow 
situation. To strengthen this statement we present in Fig. [9] the plankton 
dynamics when the upwelling is removed. We note that the concentration 
values for phytoplankton and zooplankton are slightly lower compared to 
the upwelling regime, but qualitatively there is no change observable. Thus 
localized phytoplankton blooms in vortices are possible in the wake of an 
island just due to the mechanism discussed in Subsec. [33] without any extra 
nutrient supply due to upwelling. 



4 Conclusions 

We have analyzed the interplay between hydrodynamic mesoscale struc- 
tures and biological growth in the wake of an island. Parameter values for 
the kinematic hydrodynamic flow were chosen to match the observations for 
the Canary island region, but since the basic hydrodynamic features studied 
here are commonly observed in other areas too, we expect our results to be 
of general validity. Our study is focused on the emergence of a plankton 
bloom localized in a vort ex in the wake of an island. In a previous paper 
(jSandulescu et al.l . boOTbl ) it has been pointed out that under certain con- 



ditions a vortex may act as an incubator for plankton growth and primary 
production. Here we have revealed the mechanism of such a plankton bloom. 
If the hydrodynamic flow far away from the island is dominated by a jet. 
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then the hydrodynamic time scale is much faster than the biological one, so 
that considerable growth of plankton cannot be observed. By contrast, in 
the wake of an island we obtain a much slower time scale which becomes 
comparable to the biological one giving rise to an exponential growth of phy- 
toplankton and thus to the emergence of a plankton bloom within a vortex. 
The essential factors for this phenomenon to happen are (i) the long resi- 
dence times in the vicinity of the island leading to an enrichment of nutrients 
and plankton in the neighborhood of the island; (ii) the transport and subse- 
quent entrainment of nutrients and plankton to the interior of the vortex due 
to filamental structures emerging with the chaotic saddle beyond the island, 
and (iii) the confinement of plankton in the vortex. Though the upwelling of 
nutrients in an upwelling zone enhances the emergence of localized plankton 
blooms, it is not a precondition for this phenomenon to occur. The extra 
nutrients supplied by vertical mixing in areas away from the vortex street 
are not a part of the mechanism explained here. Upwelling could be more ef- 
fective if the vortices directly cross upwelling zones whe n travelling through 



the ocean. Similar si tuations have been considered in (jMartin et al.l . |2002| : 



Pasquero et al. . 20051 ). There it has been shown that under conditions where 



upwelling occurs only directly within vortices, a plankton bloom within a 
vortex can be initiated. Therefore this mechanism, which relies mostly on 
upwelling, is different froin the o ne reported here. The variety of real ob- 
servations (jAristegui et al.l . 119971 ) in the Canary wake may benefit from the 
identification of the different possible mechanisms. 



Appendix A: The numerical algorithm 

The investigation of the interplay of biological and physical processes is 
based on advection-reaction-diffusion systems (Eqs. [6]). This system of 
partial differential equations is solved numerically by means of a semi- 
Lagrangian algorithm. The concentration fields of nutrients A^, phytoplank- 
ton P and zooplankton Z are represented on a grid of [500 x 300] points. 
The integration scheme splits the computation into three steps correspond- 
ing to advection, reaction and diffusion which are performed sequentially in 
the following way: 

1. Advection: Each point of the grid is integrated for a time step dt 
backwards in time along the trajectory of a fluid parcel in the velocity 
field. This procedure yields the position from which a fluid parcel 
would have reached the chosen grid point. Typically this position is 
not located on a grid point but somewhere in between. 
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2. Reaction: Once the position of the fluid parcel in the past is found, 
we compute the values of the concentration fields of N, P and Z at 
this point and take them as initial values for the reaction term which 
is integrated forward in time for a time step dt. Since the position 
of the fiuid parcel is not on a grid point the concentration fields have 
to be evaluated by means of a bilinear interpolation using the nearest 
neighbor grid points. 

3. Diffusion: Finally we perform a diffusion step based on an Eulerian 
scheme. Note, that the reaction step induces already a numerical 
diffusion of the order oc dx^ /dt due to the interpolation. Therefore 
one has to make sure that the real diffusion according to th e Qkubo 
estimate ( D = lOm^/s, dimensionless value D = 0.041472) (jOkubol . 



Il97lh is larger than this numerical diffusion. Additionally the stability 
condition of the Eulerian diffusion step {Ddtd/dx^ < 1 with dt^ the 
diffusion time step) has to be fulfilled. Both conditions together require 
that the diffusion time step dtd is much smaller than dt. We have 
chosen dt^ = dt/ 10 in our algorithm, i.e. after each step of advection 
and reaction we perform 10 steps of diffusion. The parameters used in 
the computation are dx = 0.02, dt = 0.01 and dtd = 0.001 expressed 
in units of = 30 days for time and r = 25 km for space. 

Appendix B: Finite Size Lyapunov- Exponents 

Stretching by advection in fluid flows is often described by means of Lya- 
punov exponents. They are defined as the average of the exponential rate of 
separation of initially infinitesimally separated parcels. For application with 
data sets from tracer experiments the infinite time limit poses a problem. 
To study non-asymptotic dispersion pr ocesses, Fini t e Siz e Lyapunov Ex- 



pone nts (FSLE) have been introduced (jArtale et al.l . 119971 : Id'Qvidio et al 



2OO4I ). The FSLE technique allows us to characterize dispersion processes 



and to detect and visualize Lagrangian structures, such as barriers and vor- 
tices. The FSLE are computed by starting two fluid elements at time t close 
to the point x but at a small distance 5q, and let them to evolve until their 
separation exceeds 6f. From the elapsed time, r+, the FSLE is calculated 
as 

A+(x,f,5o,(5/) = — log/ (7) 

The positive subindexes indicate that the tracers are advected forward in 
time. A+ is a scalar measure giving the stretching rate in the flow as it is 
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the inverse of the separation time r. 

The same definition can be applied to tracers integrated in the negative 
direction in time. A_ gives the contraction rate in the flow at the specified 
position: 

A_(x,i,5o,<^/) = — log§^ (8) 

Regions with high values of A+ and A_ trace out approximately the stable 
and unstable manifolds of the chaotic saddle. These manifolds cannot be 
crossed by fluid parcel trajectories and thus greatly influence the transport 
in the area. 
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Figure 2: The concentration of phytoplankton (normalized with the 
steady state concentration, Pamb) for inflow at ambient concentrations (left) 
and for low inflow concentrations (right). Snapshots taken at t/Tc = 
4.1,4.35,4.6,4.85. 
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Figure 3: The time evolution of the concentrations (C) of the species and 
the primary productivity. Starting concentrations are 0.01 x steady state 
concentrations of N, P and Z (top)2§nd 0.3 x steady state concentrations 
of N, P and Z (bottom). 




Figure 4: The trajectory of a fluid parcel released in the flow at the coordi- 
nates (-1.15, 0) at time t = 0. Its subsequent positions are plotted with a 
dimensionless time step At = 0.001. 
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Figure 5: The unstable manifold of the chaotic saddle. Snapshots taken at 
t = 1.5 Tc (lower plot), and 2 Tc (upper plot). 
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Figure 6: Residence time (coded in color in units of Tc) of fluid elements 
inside the area as a function of its initial positions at time t/Tc = 0.25. 
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Figure 7: The number N of tracers inside the area [0, 2] x [— 1, 1] as a function 
of time in units of the period of the flow Tc. Vortex strength w = 200. Left 
panel shows the very early decay, with an exponential fit. Right panel plots 
the overall decay, with a power law fit at large times. 
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Figure 8: Finite Size Lyapunov Exponent distributions. We plot the field 
A+ — |A_| at time t/Tc = 0.25. Stable and unstable manifolds in the flow 
are approximated by the most positive and most negative filaments in the 
distributions. 
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Figure 9: The concentration of phytoplankton in the absence of upwelhng. 
Snapshots taken at t/T^ = 4.1,4.35,4.6,4.85. The color coding is in the 
same range as in Fig. [2j 



